Photon solid phases in driven arrays of nonlinearly coupled cavities 
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We introduce and study the properties of an array of QED cavities coupled by nonlinear elements, 
in the presence of photon leakage and driven by a coherent source. The nonlinear couplings lead to 
photon hopping and to nearest-neighbor Kerr terms. By tuning the system parameters, the steady 
state of the array can exhibit a photon crystal associated to a periodic modulation of the photon 
blockade. In some cases the crystalline ordering may coexist with phase synchronisation. The class 
of cavity arrays we consider can be built with superconducting circuits of existing technology. 
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Introduction. — Since its beginning, the study of light- 
matter interaction in cavity and circuit QED has been 
providing a very fertile playground to test fundamental 
questions at the heart of quantum mechanics, together 
with the realisation of very promising implementations of 
quantum processors [TJ [2] . The coupling of separate cavi- 
ties through photon hopping introduces an additional de- 
gree of freedom that is receiving increasing interest both 
theoretically and experimentally. 

Cavity arrays, periodic arrangements of neighbouring 
QED cavities, have been introduced [3H5] as prototype 
systems to study many-body states of light. Their very 
rich phenomenology arises from the interplay between 
strong local nonlinearities and photon hopping. In the 
photon blockade regime the array enters a Mott insu- 
lating phase, where photon number fluctuations are sup- 
pressed. In the opposite regime, where the hopping dom- 
inates, photons are delocalised through the whole array 
with long-range superfluid correlations. The phase dia- 
gram has been thoroughly studied by a variety of meth- 
ods and the location of the different phases, together with 
the critical properties of the associated phase transitions, 
have been determined (see, e.g., the reviews [SHE])- 

The properties of cavity arrays resemble in several as- 
pects those of the Bose-Hubbard model [9], as long as 
particle losses can be ignored. Cavity arrays, however, 
will naturally operate under nonequilibrium conditions, 
i.e. subject to unavoidable leakage of photons which are 
pumped back into the system by an external drive. In 
this case the situation may change drastically, and, to a 
large extent, it is an unexplored territory. Only very re- 
cently the many-body nonequilibrium dynamics of cavity 
arrays started to be addressed [T0UT3] , thus entering the 
exciting field of quantum phases and phase transitions in 
driven quantum open systems [T4HT9] . 

Since the very beginning, all the works devoted to cav- 
ity arrays studied the case in which adjacent cavities are 
coupled by photon hopping. In this Letter we introduce 
a new class of arrays in which the coupling between cavi- 



ties is mediated by a nonlinear element /medium. Thanks 
to the flexibility in the design of the nonlinear coupling 
elements, these finite-range couplings can appear in the 
form of cross-Kerr nonlinearities and/or as a correlated 
photon hopping, leading to a steady-state phase diagram 
that is a lot richer than the cases which have been con- 
sidered so far. Here we discuss in particular the appear- 
ance of a new phase in cavity arrays, a photon crystal, 
which emerges in the steady-state regime when the array 
is driven by a coherent homogeneous pump. 

A technology that is very well suited for realizing cavity 
arrays with such features is provided by circuit QED [20] , 
where exceptional light-matter coupling has been demon- 
strated [21], first experiments with arrays of up to five 
cavities have been done [22], and great progress towards 
experiments with lattices of cavities has already been 
achieved [8]. 

In the following we first introduce the model for the 
cavities with their nonlinear couplings, the external drive 
and the unavoidable leakage of photons. We present a 
possible implementation in an array of circuit- QED cav- 
ities that are coupled via a nonlinear element. We then 
study the steady-state regime by means of a mean-field 
approach and Matrix Product Operator (MPO) simula- 
tions [23]. The scenario that emerges is rather complex, 
with the appearance of a number of phases and phase 
instabilities. We focus in particular on the possibility of 
spatial photon patterns that can emerge. For bipartite 
lattices the photon blockade is modulated on two dif- 
ferent sublattices, furthermore on increasing the photon 
hopping it may also coexist with a global coherent state. 

The model. — The cavity array is sketched in Fig. 
The coupling between the cavities is mediated by a non- 
linear element. In the specific implementation in circuit- 
QED this element is a Josephson nano-circuit. When 
the coupling between the cavities is realized through the 
circuit described in Fig. [TJd, linear tunneling of photons 
between adjacent cavities can be tuned and even fully 
suppressed by adjusting the nonlinear coupling circuits to 
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FIG. 1: (color online), a) An array of QED-cavities described 
by oscillator modes (red circles) that are coupled via nonlinear 
elements (crossed boxes), b) and c) Implementation of its 
building blocks in circuit-QED for one- and two-dimensional 
lattices. The circuit cavities are represented by a LC-circuit 
with capacitance C and inductance L and mutually coupled 
through a Josephson nano-circuit, with capacitance Cj and 
Josephson energy Ej, that generates the on-site and cross- 
Kerr terms in Eq. M. Details of this implementation can 
be found in Ref. [24] . An alternative approach to cross-Kerr 
interactions is discussed in Ref. 1251. 



a suitable operating point, see Supplementary Material 
(SM) [24 for details. In this regime the cavities are cou- 
pled via a strong cross-Kerr term and further correlated- 
hopping terms, which can lead to considerable modifi- 
cations in the phase diagram [25]. Yet there are also 
more involved approaches involving multiple transmon 
qubits to realize cross-Kerr interactions in the absence of 
correlated hoppings [26], see also [27l [28]. However, in 
the regime of parameters we are interested in and where 
a photon crystal emerges, correlated hopping leads to 
small quantitative corrections at the expense of compli- 
cating considerably the analysis. In the SM we quantify 
these differences in more details [24] . 

Here for simplicity, we concentrate on the salient fea- 
tures of the array in Fig. [T] which are captured by the 
effective Hamiltonian (in the rotating frame), 

U = ^[-5a\ai + n(a\+ai)]- J ^(aja, + H.c.) 

i (id) 

+ U J2 n i( n i ~ !) + V J2 nifl i ' W 



where the number operator rii = a\ai counts the pho- 
tons in the z-th cavity (aj/a^ being the corresponding 
creation/annihilation operators). The first three terms 
describe respectively the detuning 5 of the cavity mode 
with respect to the frequency of the pump, the coherent 



pump with amplitude ft and the hopping of photons be- 
tween neighboring cavities at rate J. The last two terms 
take into account the nonlinearities through the onsite- 
and cross-Kerr terms with the associated energy scales U 
and V respectively. In the specific case of circuit-QED 
arrays, the two types of nonlinearities can be realized 
through the setup of Fig. [I] In order to keep our results 
as general as possible, we consider the effective model ([I]) 
without specifying further the underlying matter-light in- 
teraction term. 

The dynamics of the array is governed by the Master 
equation 



p = -i[H, p] + | ^2(2(npal 



nip - prii 



(2) 



where ft -1 is the photon lifetime in each cavity. The 
model in Eq. ([!]) together with Eq. Q encompasses, in 
some limiting cases, regimes that were already addressed 
in the literature. The regime of U — » oo and J = was 
considered in Ref. [18 , where an antiferromagnetic phase 
was first predicted in Rydberg atoms. The case of on-site 
Kerr nonlinearity, i.e. V = 0, is the only one studied so 
far in cavity arrays [12] . The model considered here offers 
a much richer phase diagram. A unique characteristics 
of the cavity arrays with nonlinear couplings is that the 
cross-Kerr nonlinearity V can even exceed U . By cou- 
pling an additional qubit locally to each resonator [29] , 
different ranges of the ratio V/U can be explored. More- 
over, in devices where on-chip control lines can be used 
to locally thread magnetic fields through the loops of 
the coupling circuits and the additional transmons, the 
ratios J/U respectively J/V and V/U can be tuned on- 
chip. For this reason we will, in the following, consider 
J7, V and J as independent. 

We first discuss the steady-state phase diagram in the 
mean-field approximation, which becomes accurate in 
the limit of arrays with large coordination number z. 
The decoupling in Eq. is performed on the hopping 
and cross-Kerr terms, z~ x a l a j ~> ( a \) Y^ieB a i + 

(4)E iG A a i and £<t,i> ( n A)T,ieB n i + 

( n B)Yj e A n 3i wnere we assumed a bipartite lattice, 
A and B being the two sublattices. The mean- field 
analysis simplifies the dynamics dictated by Eq. ^ 
to two coupled equations for the two different sublat- 
tices. As a function of all the parameters characterizing 
the system and its dynamics, one gets a very rich be- 
havior in the asymptotic regime which includes steady- 
state/oscillating phases, as well as uniform/ staggered 
configurations. Here we highlight what we think are its 
most intriguing features. All the couplings will be ex- 
pressed in units of the photon lifetime, k = 1. 

Mean-field steady-state diagram. — Fig. [2j where for 
the moment we set the hopping to zero, shows that, 
on increasing the cross-Kerr term, the array can reach 
a steady state in which the photon number is modu- 
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FIG. 2: (color online). Order parameter An for the photon 
crystal in the U — V plane at zero hopping. If the cross- 
Kerr term exceeds a critical threshold V c , the steady state is 
characterized by a staggered order in which An 7^ 0. Here we 
fixed Q = 0.75 and S = 0, for which zV c « 0.44 at U = 0, 
while zV c ~ 5.73 in the hard-core limit (U —> 00). In the 
inset we show An as a function of Q and V at a fixed value 
of U = 1. Here and in the next figure the color code signals 
the intensity of the order parameter, while dashed green lines 
are guides to the eye to locate the phase boundaries. 

lated as in a photon crystal, the order parameter be- 
ing An = \(tia) — (^b)I- Here the area above the 
green line denotes the crystalline phase. In the U — >• 00 
limit, the transition between the uniform and the crys- 
tal phase is located at zV c ~ 7^ (—2 S + • v /7 00 )/4^ 2 with 
7oo = 4 J 2 + SQ 2 + 1, the previous expression holding 
for small detunings, and coincides with the transition to 
an antiferromagnetic phase described in Ref. [T8] . Note 
that a lower value of U favors the crystal phase. In the 
opposite limiting case of U = 0, the transition is found 
at zV c ~ 7o(-2<5 + y^)/^ 2 with 7o = 4£ 2 + 1. We 
deliberately considered a regime in the parameter space 
where V > U since, as already mentioned, it is a peculiar 
feature of the cavity arrays proposed here. The transi- 
tion to the crystal phase is reentrant as a function of the 
drive (inset to Fig. [2J. At very small pumping the den- 
sity is too low to lead to a photon crystal. Vice-versa it 
also disappears on increasing since pumping favors an 
homogeneous photon arrangement. A similar feature has 
been observed in the limit U = 00 [T8] . 

If the hopping between photons is switched on, der- 
ealization will suppress the solid phase and at a critical 
value of J (which depends on V, U, 5, and Q) there is a 
transition to a normal phase. This is shown in Fig. [3J as 
a function of the cross-Kerr nonlinearity V (Fig. [3^t) and 
of the detuning 5 (Fig. [3)3). In Fig. |3^i we display the 
case [7=1, while at smaller values of U the phase dia- 
gram shows a reentrance. Although interesting, further 
analysis is needed to see if this feature is only present at 
mean-field level. Yet it does not seem improbable that 
an increased hopping could facilitate the redistribution of 
particles into a crystalline order imposed by the interac- 
tions. We conclude this discussion by pointing out that, 
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FIG. 3: (color online). Order parameter for the photon crys- 
tal at finite hopping. An is plotted as a function of J and 
V, for 5 = and U — 1 (panel a), and as a function of J 
and 5, for zV = 0.6 and U — (panel b). Here we fixed 
^1 = 0.75. At finite values of the detuning, in addition to the 
normal (white) and crystalline (colored) phases, an interme- 
diate region (shaded green), characterized by an oscillatory 
behavior in the asymptotic state, appears. As discussed in 
the main text, we suggest this last regime may be seen as a 
nonequilibrium analog of a supersolid. 



as discussed in Ref. [24 , under nonequilibrium condi- 
tions it is even possible, although much harder, to realize 
a crystalline phase at V = 0. As a matter of fact in that 
case, for some values of the coupling constants, the steady 
state can be either uniform or crystalline, depending on 
the initial conditions. 

The mean-field phase diagram in the 5 — J plane and 
for U = is depicted in Fig. (3}d. As highlighted in the 
region between the dashed green lines, for 0.8 < S < 1.1, 
on switching on the photon hopping, a new intermediate 
phase appears. In this region, even in the long-time limit 
the state never becomes completely stationary and there 
is a residual time-dependence of (a) with (a a) 7^ (ob)? 
i.e. there is an additional time dependence of (a) on 
top of the trivial oscillation with the frequency of the co- 
herent drive that is hidden in our choice of the rotating 
frame. At the same time the system shows An ^0. In 
Fig. [4^i we show the time evolution of the real and imagi- 
nary parts of (a) for the two sublattices. A closer inspec- 
tion of the properties of the oscillating phase reveals that 
the reduced density matrix of a single site (in either of 
the two sublattices) is a coherent state which evolves pe- 
riodically in time as shown in Fig. ^p. There we plotted 
the Wigner function W(x,p) = J(x — jj\pa\x + y)e 2xpy dy 
of one sublattice at a given time, with x = (a + af)/y/2, 
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FIG. 4: (color online), a) Time dependent traces of the real 
and imaginary parts of (a) for the two sublattices as a function 
of time in the steady-state regime of the intermediate oscillat- 
ing phase of Fig. [3)3. b) The Wigner transform of the reduced 
single-site density matrix for sublattice A is plotted at a given 
time in the intermediate oscillating phase of Fig.|3]3. Here we 
used the same parameters as in Fig. [3]o and fixed z J = 0.2, 
£ = 0.9. 
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p = i(o) — a)/y/2, and \x) being an eigenstate of the po- 
sition operator x. Following the analysis performed in 
Ref. [19 j, we are lead to conclude that in this region the 
dynamical evolution of the whole array is synchronised 
separately in the two different sublattices. The contem- 
porary presence of checkerboard ordering and global dy- 
namical phase coherence suggests us to view this inter- 
mediate phase as a nonequilibrium supersolid phase [30 . 
The intermediate region extends also at finite- U values, 
although the coherent state of Fig. will be progres- 
sively deformed on increasing the on-site repulsion. 

MPO simulations. — Most of the features we discussed 
can already be seen in small arrays. To show some ex- 
amples and to further support the mean-field analysis 
given above, we here present results that were obtained 
for linear chains of cavities, z = 2, with MPO simu- 
lations [12] of the Master equation ([2|, which provide 
a (numerically) exact description of its nonequilibrium 
many-body dynamics. Fig. [5^i shows the density-density 
correlation function g^ 2 \i,j) = {a\a)-ajai) / (ni){rij) for a 
chain of 20 cavities with S = 0, J = 0, U = 0.5, Q = 0.4 
and various values of V. One clearly sees that a stag- 
gered dependence of the distance r = \i — j|, indicating 
strong density-density correlations, appears for nonzero 
V, whereas for V = photons in distinct cavities are un- 
corrected (</ 2 ) = 1). Fig. JHJd shows for a chain of 21 
cavities with £ = 0, U = 1, V = 1, O = 0.4 and various 
values of J. The spatial range of density-density correla- 
tions shrinks with increasing tunneling rate J, indicating 
a crossover to an uncorrelated state. A more quantita- 
tive analysis of the decay of correlations with increasing 
distance is not conclusive for the chain length considered 
here. A true ordering in the steady state can probably 
only be stabilized in two-dimensions. 

Conclusions. — In this Letter we introduced cavity ar- 
rays with coupling mediated by nonlinear elements. This 



FIG. 5: (color online). MPO results for linear chains of cav- 
ities, a) g (2) (10J) for 8 = 0, J = 0, U = 0.5, Q = 0.4 and 
V as in the legend, b) g (2) (llj) for 5 = 0, U = 1, V = 1, 
Q — 0.4 and J as in the legend. 



opens the way to study a variety of new possibilities, in- 
cluding correlated photon hopping and finite-range pho- 
ton blockade. We concentrated on this last point study- 
ing the effect of a cross-Kerr nonlinearity on the steady 
state and found a very rich phase diagram. A photon 
solid characterized by a checkerboard ordering of the av- 
erage photon number appears for a substantial range of 
the coupling constants. In addition we see that, for some 
choice of the parameters, a finite hopping stabilizes a 
phase where the crystalline ordering coexists with a glob- 
ally synchronized dynamics of the cavities, suggesting an 
analogy to a nonequilibrium supersolid. Most of the re- 
sults presented in this work were obtained in a mean- 
field approximation. We corroborated the existence of a 
steady-state solid phase by studying a one-dimensional 
array by means of a Matrix Product Operator approach. 
This last analysis confirms that a crystalline ordering of 
photons can be observed with existing experimental tech- 
nology. 
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Supplementary Material 



Cross-Kerr nonlinearities in arrays of circuit cavities 



Here we describe a circuit quantum electrodynamics (cicuit-QED) setup that represents an array of cavities coupled 
by nonlinear elements and can be effectively described by a lattice of harmonic oscillators coupled via cross-Kerr 
nonlinearities. The setup thus provides an experimentally feasible way for implementing the model described by 
equations (1) and (2) of the main text. Cicuit-QED is particularly well suited for this task because of the great design 
flexibility, the tunable nonlinearity provided by Josephson junctions and the exceptionally high coupling between 
subsystems that can be reached. 

For our derivation, we consider lumped element resonators, representing the cavities in our array, which are coupled 
conductively through capacitively shunted Josephson junctions, see Fig. [6] for a sketch of the circuit representing a 
linear chain of cavities. We focus on lumped element resonators, in order to keep the derivation simple and transparent. 
Coplanar waveguide resonators work equally well. For the considered combination of capacitive and inductive coupling 
between the resonators, the linear parts of the couplings can cancel each other for suitable choices of the parameters, 
leaving a residual coupling via the nonlinearity of the Josephson junctions. 

In terms of the node fluxes fa, the Lagrangian of the whole setup reads, 



C 



E 



c 



2L 



E 

(id) 



j)j) 2 + Ej cos 



(3) 



with L and C the inductance and capacitance of the lumped element resonators and Cj and Ej the capacitance 
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and Josephson energy of the Josephson junctions. <j>o = h/(2e) is the reduced flux quantum. Assuming Cj/C <C 1 
we invert the capacitance matrix to first order in Cj/C for performing the Legendre transformation to obtain the 
Hamiltonian 



H 



nl 




(4) 



where C = C + 2Cj, l/L = 1/(2L) + l/Lj, and Lj = <f%/Ej. 

The charges on the islands qi = C<pi — Cj + O'(^) denotes all connections to site z), defined by the 

coupling capacitances of the Josephson junctions and the lumped element resonators, and the fluxes <j)i associated 
to the phase drop at the inductance of the lumped element resonators are our canonically conjugate variables. We 
quantize the Hamiltonian by introducing bosonic lowering and raising operators ai and a\ according to, 




The quantized coupling Hamilton operator thus reads, 



(5) 
(6) 



H c = y^ j ujX-(a\aj + a^aj) + ujX + (a\a\ +a^aj), 



2L 



-. We choose Cj/(C+2Cj) = 2L/(2L+Lj) 



where the oscillator frequency is uj = 1/ V LC and X± — 2 l+Lj c+2Cj 
such that X_ = and the linear tunneling between neighboring oscillator sites vanishes. Furthermore we neglect the 
terms proportional to X + in a rotating wave approximation which is justified as long as Cj/(C+2Cj)+2L/ (2L+Lj) <C 
2 holds. In this way we arrive at a chain of harmonic oscillators that are decoupled in linear order. 

The nonlinear parts of the Hamilton operator however provide us with coupling terms between the neighboring 
oscillators and on-site nonlinearities. We restrict ourselves to fourth order nonlinear it ies, perform a rotating wave 
approximation and arrive at the many-body Hamiltonian, 



H 



E 



(uj + 5(jj)a\cii — az ~^~ 



(xz—a]a]ajaj 



(id) 



rr(ch) 
H i,3 



(7) 




Site i 



FIG. 6: (color online). Electrical circuit sketch of the setup we envision to realize a system with cross Kerr nonlinearities. 
Here for the example of a one-dimensional chain of cavities. Lumped element LC-circuits are coupled via capacitively shunted 
Josephson junctions. The Josephson energies could be tuned in situ by replacing the Josephson junctions with superconducting 
interference devices. 
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where Sou is a small correction coming from the normal ordering process of the nonlinearity and Eq — e 2 /(2C 2 ) is 
the charging energy of the individual lattice sites of our model. The second term in Eq. ^ is an on-site Kerr term 
whereas the third term describes cross-Kerr interactions. This specific implementation thus realizes the Hamiltonian 
of equation (1) in the main text for U = —azEc/2 and V = — 2aEc, thus leading to zV/U = 4. 

The term h[°^ = aEcfaiOjCijCij + a\a\aiCLj — ^ a \ a \ a j a j) + H.c. describes correlated hopping of photons between 
adjacent sites i and j. For the purposes of the present work, i.e. a description of solid and supersolid phases, this 
term can be treated at a mean- field level. A mean- field approximation would replace a\ai —> (aja$) and the first 
two terms will thus only give rise to a (possibly sublattice dependent) renormalization of the single-photon hopping. 
The term ajaja^ + ia^ + i in turn is a factor of 4 weaker than the cross-Kerr term. Moreover, in regimes of large Kerr 
interactions, double occupations of lattice sites are small and this term becomes ineffective. For these reasons, in 
most of our analysis we disregard the correlated hopping at present. A signature of the effects induced by H- ^ is 
evident from the data shown in Fig. [7| where we plotted the phase boundary between normal and crystalline phase 
in the U — V plane at zero hopping. The relatively small discrepancies in the curves with and without such terms 
points towards the fact that only quantitative modifications are induced in the parameter space we are considering, 
while all the qualitative features of our results should be unaffected. It should however be noted that a more accurate 
treatment may lead to new phases in the diagram such as a pair supernuid state, when the corresponding hopping 
term becomes comparable with the cross-Kerr nonlinearity. This is however beyond the scope of the present analysis 
and it will be treated separately. 
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FIG. 7: (color online). Emergence of the photon crystal in the U — V diagram at J = 0. We highlight the effects induced 
by the terms in by comparing mean-field results without (green line) and with (blue line) such contributions. Here the 

strength of the correlated-hopping terms is set equal to —U, while we set Q = 0.75 and S = 0, as in Fig. [2] of the main text. 

Our approximations require that a <C 1. Nonetheless the cross-Kerr interaction can be much larger than photon 
losses, i.e. 2aEc ^> ft, since, e.g., transmon qubits have Ec/h ~ 0.5GHz and T\ ~ 1/is [2|. 

Finally let us point out that the Josephson junctions linking two neighboring oscillators can be built tunable by 
replacing them with a dc- SQUID. In this way the Ej and thus the Lj can be modulated by applying an external flux 
to the dc-SQUIDs and the Hamiltonian ([7]) can be tuned in real time. Hence, by choosing the external flux such that 
X_ / 0a linear tunneling of photons between the resonators can be introduced. Moreover larger on-site nonlinearities 
can be introduced by coupling each resonator locally to a superconducting qubit, e.g. a transmon. 

Photon crystal at V — 

We conclude the discussion of the solid phase by showing that the realization of a steady-state photon solid does 
not necessarily requires a cross-Kerr nonlinearity. Surely the presence of V ■=£ stabilizes the solid in a wide range of 
the parameter space, thus making it easier to be observed experimentally. Nevertheless, while in equilibrium V ^ 
is a necessary requirement, under nonequilibrium condition this ceases to be the case. The only requirement is a 
finite photon hopping and an initial unbalance in the occupation of the two sublattices. In Fig. [5] the regions in which 
the steady-state solid phase is reached are plotted as a function of the initial unbalance. Note that the value of the 
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FIG. 8: (color online). Phase diagram at V — in the space of initial conditions, with zJ — 6.2, 8 — 0.8, Q = 2 and in the 
hard-core limit U — ► oo. White areas denote starting occupation values where a steady-state uniform phase is reached, while 
red areas denote initial conditions leading to an antiferromagnetic phase with An ~ 0.1027. 



crystalline order parameter An does not depend on the choice of the initial conditions. 
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